Mutual Inductance Route to Paramagnetic Meissner Effect in 2D Josephson- Junction 

Arrays. 
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We simulate two-dimensional Josephson junction arrays, including full mutual-inductance effects, as 
they are cooled below the transition temperature in a magnetic field. We show numerical simulations 
of the array magnetization as a function of position, as detected by a scanning SQUID which is placed 
at a fixed height above the array. The calculated magnetization images show striking agreement 
with the experimental images obtained by A. Nielsen et ala The average array magnetization is 
found to be paramagnetic for many values of the applied field, confirming that paramagnetism can 
arise from magnetic screening in multiply-connected superconductors without the presence of d- wave 
superconductivity. 
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A DC paramagnetic susceptibility, reported first by 
Braunisch et alB for BSCCO, occurs in many high — T c 
superconductors when cooled through their transition 
temperature in an external magnetic field. This sur- 
prising result, known as the paramagnetic Meissner ef- 
fect (PME) , contrasts with the standard diamagnetic re- 
sponse of classical superconductors and has been subject 
of extensive investigations for the last ten years. 

Some theoretical workEl suggested that the PME pro- 
vided indirect evidence for d-wave symmetry in the super- 
conducting order parameter. In this picture, 7r-junctions 
formed between misaligned grains were the cause of the 
anomalous magnetic response. 

PME observed in low — T c superconductors with s- 
wave order parametersu demonstrated that 7r-junctions 
were not required for PME. New theories for PME were 
developed, advocatiag non-equilibrium phenomena such 
as flux ccanpression,0 surface barriersu and a giant vor- 
tex state.0 However, in the case of high — T c samples 
like BSCCO, experimental showed clearly that the gran- 
ular nature of the samples was a crucial ingredient for 
the occurrence of the phenomenon. This suggested us- 
ing arrays of (non-7r) Josephson junctions!^ as a model 
system for studying PME in granular high — T c sam- 
ples, to test whether 7r-junctions were also an essential 
ingredient. Numerical simulations of simplified Joseph- 
son junction networks (a single_multi-junction loopB or 
multi-junction concentric loopsE3) indeed showed a para- 
magnetic response. Experiments also gave indirect evi- 
dence for PME in the AC susceptibility of arrays.Eil 

Because of the many theories predicting PME in both 
s and d-wave superconductors, more stringent and de- 
tailed experimental tests were needed to find the end of 



this maze. Experiments using scanning SQUIDs were 
thus performed on high — T c superconductors!! 2 ] and on 
arrays of noLtJr junctions.!!! A scanning SQUID micro- 
scope (SSM)EJ measures the spatial distribution of the 
magnetization. The complexity of the results and the 
experimental technique pose new theoretical challenges 
in the qualitative and quantitative interpretation of the 
magnetic images. 

Here we show that a model of 2D arrays with full mu- 
tual inductance interactions captures the essential facts 
about PME in Josephson junction arrays. 
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FIG. 1. Sketch of array design. Niobium crosses are in 
two layers, light and dark grey. The Josephson junctions are 
formed at the cross overlaps, as indicated, and the external 
flux is applied perpendicular to the array. 

The arrays measured in Nielsen et al. had unit a cell 
size of 46 /im and were cooled in external flux from 
zero up to 12 $0 per unit cell of the array. A sketch 
of the array is shown in Fig. 1. The junctions had a 
J c = 600 A/cm 2 with a junction area of 5 x 5/im 2 and 
a calculated self-inductance of V — 64 pH, yielding a 
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f3 L = 2irL'I (T)/§ = 30 at 4.2 K. The experiment in- 
volved cooling the array in an externally applied field and 
then measuring the magnetization with the field still ap- 
plied. These parameters are similar to those in BSCCO 
which exhibits PMEE3 and are the parameters-used here. 

We simulate a network of N r x N c j unctions .EfTtll Using 
a vector notationJIS the current in each junction can be 
modeled by the resistively-capacitively shunted junction 
(RCSJ) model as: 



P = Iq sin (p ■ 
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Here Iq sin ip represents the current through the Joseph- 
son element (sin ip is the vector given by applying sin to 
the components of 0), and (<I>o/27t) <p is the voltage drop 
across the quasiparticle conductance G. Finally C and Iq 
are the junction capacitance and the Josephson critical 
current. 

To satisfy the Kirchhoff 's law for the currents in each 
node, we define loop currents I s connected to the junc- 
tion currents by the relationship I b = KI S (for a discus- 
sion of the role of loop currents cf. Refs. |lq , |l7|) where the 
matrix K depends on the array geometry. The fluxoid 
quantization rule for each elementary loop in the array 
gives another set of equations: 



2nP -> 
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where M performs the (oriented) sum of the phases 
around a loop; the vector / represents the normalized 
flux <f> ext — /$o due to an external field in each loop, 
i.e., the so-called frustration; ft is a vector of "quantum 
numbers" for the flux quanta in each loop; and the last 
term is the field induced by the currents flowing in all 
other loops of the array (0 mduced = L'LI S ). The matrix 
L, the mutual inductance matrix of the array (normal- 
ized to the self inductance of the single loop), represents 
the mutual coupling between loops in the arrays. Here 
we compute L by a thin wire approximation except for 
the self-inductance of a single loop(cf. Ref. Insert- 
ing the fluxoid quantization in I b = KI S , using Eq.(^j) 
we obtain a system of equations in normalized units, con- 
taining only the phase variables: 
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Here time is normalized to a cell frequency (u>~ 2 = L'C) 
and the usual Stewart-McCumber parameter appears, 
j3 c = 2irI (T)C/<boG 2 . The term m represents the nor- 
malized loop magnetization (cf. Eq. (|)). An explicit 
form for magnetization can be written as follows by in- 
verting the static form of Eq. (0) : 



m = ^-L (k T Kj 1 K T sin <p, 



(4) 



which generalizes the single-loop Eq.(l) by Nielsen et al. 
In the case of a single loop, for large (3l, there are at 
least four states which are non-degenerate and that are 
either diamagnetic or paramagnetic. The lowest energy 
states are diamagnetic for^</<£+l/2 with I integer, 
and paramagnetic for i + 1/2 < / < £ + 1. For a single 
loop, half the states are diamagnetic and half are para- 
magnetic. This contrasts with the experiments on large 
arrays by Nielsen et al. that show a clear prevalence of 
paramagnetism for / > 3. In other words, the single- 
loop model cannot explain the experimental results, even 
qualitatively. 
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FIG. 2. Simulated field cooled 10 by 40 Josephson junction 
array for a frustration / = 1.2. Parameters of simulations are 
/3l(4.2 K) = 30, /3c (4.2 K) = 66. a) Image of the array 
magnetization at z — 0; b) simulated SSM image of array 
magnetization at z = 1, sampled at positions corresponding 
to the center of array loops. The light-gray loops are the dia- 
magnetic ones, c) Histogram of loop magnetization at z = 
; d) histogram of magnetization as read by SSM at z = 1 . 

We can do a mean-field type of treatment of the tem- 
perature dependence by using the fact that (3c and Pl 
are the only temperature dependent quantities in these 
equations. Thus, we simulated field cooling in the arrays 
solving Eq. (||) for the phases and calculating the result- 
ing currents and magnetization. The simulation starts 
with a zero screening term in the equation, (3l = 0, and 
f3c = 0, as representing T > T c . Non-zero frustration / 
was fixed in the beginning of the simulation. Then, @l 
and Pc ar e increased in steps, until they reach their final, 
low-temperature value. The dynamical terms, i.e., (p and 
tp go to zero after a transient. A variable transient time 
permits control of the speed of the simulated field cooling 
process. We used parameters similar to the experiments,!!] 
i.e., Pl{T = 4.2 K) = 30, C (T = 4.2 K) = 66. The tran- 
sient time for each step increase in /3l ranges from 80 
to 400 normalized time units, and a typical run takes 
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30 steps. The initial conditions for the array are chosen 
with all the phases being zero and a random distribution 
of "quantum numbers" ft, simulating the disorder due 
to the initial diffusion of flux quanta, when the Joseph- 
son energy barriers are small. Details oLthe integration 
routine are described in Filatrella et alSJl 




FIG. 3. The same simulated field cooled array of Fig. 2 for 
/ = 4.8. a) Image of the array magnetization at z — 0; b) sim- 
ulated SSM image of array magnetization at z = 1, sampled 
at positions corresponding to the center of array loops. The 
light-gray loops are the diamagnetic ones, c) Histogram of 
loop magnetization at z — ; d) histogram of magnetization 
values as read by SSM at z = 1. 

In order to have a significant comparison between the 
numerical simulations and the experiments, we take into 
account the SQUID-sample separation at a non-zero dis- 
tance z above the array. Typical values of z have been 
chosen within the limits indicated by Ref. |l|, 40 to 60 jtim, 
and we normalized z to the array unit cell size, 46 /um. 
The field at a distance z was built by superposition of 
the fields generated by the currents. Each current injthc 
array is modeled using the thin-wire approximation.^ 

Next, the flux within a square corresponding to the 
SQUID area was calculated, for different positions above 
the array. We chose to calculate positions corresponding 
to the centers of the array loops (i.e. one point per loop) 
at distance z above them. 

Fig. 2 reports the field-cooled magnetization for a 
10 x 40 array with / = 1.2 and clearly shows a diamag- 
netic behavior both locally and in the average magneti- 
zation. Figs. 2a and 2b respectively show the magneti- 
zation at z — and z = 1. Figs. 3 and 4 show the same 
array for / = 4.8 and / = 12.2: Figs. 3a and 4a report 
the magnetization at z = 0, Figs. 3b and 4b magnetiza- 
tion at z = 1. For values of frustration above 3, the array 
shows an overall paramagnetic response. It is interesting 



to note that in all cases, at z = 0, there is a connection 
between the simulated arrays and the simple single loop 
picture. If, for a given value of frustration, an isolated 
loop is diamagnetic (lowest energy state), for the same 
value of frustration the simulated array shows a larger 
number of diamagnetic loops. These diamagnetic loops 
form a "sea" in which a few paramagnetic loops stand 
out (cf. Figs. 2a, 4a). If the isolated loop is paramag- 
netic, the "sea" is formed by paramagnetic loops with 
few diamagnetic loops in the array (cf. Fig. 3a). 
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FIG. 4. The same simulated field cooled array of Fig. 2 for 
/ = 12.2. a) Image of the array magnetization at z = 0; b) 
simulated SSM image of array magnetization at z = 1, sam- 
pled at positions corresponding to the center of array loops. 
The light-gray meshes are the diamagnetic ones, c) Histogram 
of loop magnetization at z — ; d) histogram of magnetiza- 
tion values as read by SSM at z = 1. 

At z — 1 the mixing of flux lines produces a smeared 
flux distribution that is very similar to the experiments 
(cf. Ref. |l|) . We note that for large frustration values (cf. 
Fig.s 4a and 4b), due to different magnetization strength, 
the far-field array image is paramagnetic, although the 
corresponding state for an isolated loop is diamagnetic. 

In Figs. 2c, 3c and 4c, histograms of the loop magne- 
tization are reported at 2 = 0. We find two peaks repre- 
senting the diamagnetic ($tot — Q e xt < 0) and paramag- 
netic (<&tot — §ext > 0) loops. The peak position essen- 
tially corresponds to single loop values for the same frus- 
tration. The peak width is determined by mutual induc- 
tance effects. Generally only two magnetization peaks 
are found, one diamagnetic and one paramagnetic (with 
the exception of a few loops in the / = 4.8 case, which 
show a higher value of paramagnetic magnetization, cf. 
Fig. 3). The majority loops magnetize weakly whereas 
the minority loops magnetize strongly. Figs. 2d, 3d, and 
4d show the histograms evaluated at z = 1. Similarly 
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to the measured images, we observe a smearing effect: 
Histogram peaks merge, so that the overall distributions 
appear similar to the experimental ones. Merging of his- 
togram peaks starts approximatively at z ~ 0.3. The 
results discussed for the Fig^. 2, 3 and 4 can be extended 
to other frustration values:E3 Simulations show that for 
£</<£ + 1/2, with £ integer, the diamagnetic loops 
predominate in number, whereas lav £+1/2 < / < £+ 1 
the paramagnetic loops dominate. For / = £+1/2 the so- 
lution tends to have an equal number of diamagnetic and 
paramagnetic loops. The magnetization strength shows 
a more subtle behavior: For £ < f < £+1/2 the strongest 
magnetization is paramagnetic, for £+1/2 < f < £+1 the 
strongest magnetization is diamagnetic. If the frustration 
equals a half integer, / — £ + 1/2 the paramagnetic and 
diamagnetic peaks are of equal strength, so their average 
magnetization measures zero. 
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FIG. 5. Dependence of mean array magnetization on frus- 
tration for a 10 by 40 array. Parameters of simulations are 
/3 L (4.2 K) = 30, /3o(4.2 K) = 66, z = 1. 

In Fig. 5 we report the mean magnetization over a 
10 x 40 array, for different frustrations, at z = 1. The 
mean magnetization depends on the blend of paramag- 
netic/diamagnetic strength in the loops and their num- 
ber. A trend shifts the array magnetization toward para- 
magnetism, starting from / > 3. The mean magnetiza- 
tion depends weakly on noise: A test with different ran- 
dom distributions of quantum numbers shows that this 
accounts for an error of about 1.2%. Our estimation of 
z adds another source of error, but our simulations show 
that this error accounts for no more than 5%, with a z 
variance of 20%. On the other hand, the mean magneti- 
zation depends on the dimension of the array. A direct 
quantitative comparison with the experiments shows a 
calculated value of magnetization typically lower. Mag- 
netization strongly depends on array dimension, so the 
results presented in Fig. 5 can be only qualitatively com- 
pared with experiments, in which arrays are larger. We 
report only positive frustration (/ > 0) because Eq.(||) 
is symmetric changing the sign of frustration (the same 
array viewed from below simply maintains the same para 



and diamagnetic loops). 

We note that in all cases, i.e., both dia- and para- 
magnetic, diamagnetic behavior prevails near the array 
edges. This agrees with the experiments, which show a 
similar behavior. According to Ref. [jjthis occurs because 
the array screens the field by generating diamagnetic cur- 
rents on the array boundary and, as consequence, induces 
paramagnetic currents in the interior of the array, thus 
generating an overall paramagnetic offset. 

To further support this view, we calculated the den- 
sities of paramagnetic loops at the boundary and in the 
bulk of the array. We find that there is a clear diver- 
gence between two data sets with an increase of bulk 
density, pk with respect to boundary density, p b , for frus- 
tration / > 3. For example, at / = 1.2 the two densi- 
ties are roughly equal p b = N para /N boundary ~ 0.156 and 
Pk = N para /Nt ota i ~ 0.162, but at / = 12.2 at the bound- 
ary we have pb — 0.11 and in the bulk pk — 0.26. Tests 
on smaller arrays show that paramagnetic behavior for 
m < / < m + 1/2 arise about for N ~ 5, this is roughly 
the value predicted from Eq. (4) of Ref. [j] for (3^ = 30. 

In conclusion, the PME in Josephson junction arrays 
can be reproduced via numerical simulations which in- 
clude the full inductance matrix. The simulation re- 
sults compare favorably to experimental results: Para- 
magnetism dominates field cooling for large arrays. Sim- 
ulations also show that the single loop model is the basic 
building block describing the field cooled array behavior. 
Mutual inductance interactions create the actual distri- 
bution of loop magnetization in the arrays. The result- 
ing mean magnetization is the product of both single loop 
states and their occupancy. The observed dominant para- 
magnetism, in both experiments and simulations, arises 
from an energetic preference for paramagnetic loops in- 
terior to the array. 

Beyond this study, a number of open problems still 
remain to be analyzed. Among these, are simulations of 
larger arrays in order to make more detailed comparison 
with experiments and the study of the effect of cooling 
time and transient dynamics of the array. 
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